/* bigm_s.c - Use the solve callbacks for changing big Ms to solve the problem */

/* To run this example, command line arguments are required.
   i.e.,   bigm_s   filename
   where 
       filename is the name of the file, with .mps, .lp, or .sav extension
   Example:
       bigm_s  data/uflp_s_10x100.lp

   Make sure that the size of the problems are compatible

   Record of revisions:
   04/  2011: 
      - Output nodelps into .lp files at specified nodes
      - Show that local and global bounds can not be changed

   04/10 ~ 04/16 : 
      - Show that the bounds in the nodelp are not consistent with local bounds.
        But solutions satisifies the bounds in the nodelp rather than local bounds.
      - Make the checking routine to work and get best solutions???

   04/17:
      - Modulize by writing wrapper functions staring with DZ
      - Add logging routines for better recordings
      - The file works for problem of small size of 10x100
      - Add while loop to do check-resolve procedure
      - Introduce matrices to store the location and values of M outside of callback functions

   04/19:
      - Adding CPXgetbase routine to determine the binding constraints 
        to avoid numerical issue.
      - Set solving nodelp by simplex algorithm to garantee exact solutions
      - Explore problem of size 30x300 with startM of 1,0.99, 0.9 to obtain 
        suboptimal solutions. It was shown that optimal y are achieved but not optimal.
	Final x outputed are shown to be aborted early than necessary. 
	A guess is that the abortion is caused by fathoming or cutting off the binding nodes
	and then closing the optimality gap(Making no sense). I don't know why for the moment.
	But come back after Wednesday.
      - Show that the bounds of mip is not changed through the while loop
   04/20
      - Change CPX_PARAM_VARSEL to be 1 to avoid automatic selection by Cplex(
        which cause troubles to stop the algorithm early
      - Solidify the bigm checking routine by specify three passing criterion
          1. M_ij == 1 (By setting tolerance for difference between M_ij and 1)
	  2. x_ij == 0 (By setting tolerance for difference between x_ij and 0)
	  3. x_ij < M_ij * y_j is binding(By determing the slack variable is basic by CPXgetbase)
      - Otherwise we will change the bigm values
      - 10x100, 20x200 cases with startM of several values are tried successful
      - 30x300 case are just too slow to find solutions
      - heuristics will influece the quality of the solution, I don't why

    04/22
      - Set output to log files with size dependent names
      - Now read parameters from files instead of naively from command line
      - Single/Multiple runs can be done through bigm_single/bigm_batch shell script
      - Fix the bug of big M checking routine: decide x = 0 
        in terms of nonbasic variables
      - Make sure it works for 15x150 20x200 case.
      - Initial heuristic plays an important role for larger size case
      - (Todo) Output needs to be organized
      - (Todo) Organizde the folder
 */

#include <ilcplex/cplex.h>
//#include "cplex.h"
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>

//Add some mutex type to lock the region in the callback function
#include <pthread.h>    

/* Include declarations for functions in this program */

struct Minfo {
  int *flagPtr;
  int ni;  
  int nj;
  double stepsize;
  int **bigmRowIdx;
  int **bigmColIdx;
  double *myBestObjValPtr;
  double ***bigmValsPtr;
  pthread_mutex_t mtx;
  pthread_mutex_t mtx1;
};

typedef struct Minfo  MINFO, *MINFOptr;

static int CPXPUBLIC
   usersetM       (CPXCENVptr env, void *cbdata, int wherefrom,
                   void *cbhandle, int *useraction_p); 
static int CPXPUBLIC
  flagcallback    (CPXCENVptr env, void *cbdata, int wherefrom,
                   void *cbhandle); 

static int 
  DZturnoffpreproc(CPXENVptr env),
  DZturnoffcuts   (CPXENVptr env),
  DZsetmipstrategy(CPXENVptr env),
  DZsetmiplimits  (CPXENVptr env),
  DZsetmipoutput  (CPXENVptr env, char* logbuffer),
  DZgetmipprobsize(CPXENVptr env, CPXLPptr lp, 
		   int *nRowsPtr, int *nColsPtr, int *nBinsPtr),
  DZgetbigmindex   (CPXENVptr env, CPXLPptr lp, int ni, int nj, 
		    int **bigmRowIdx, int **bigmColIdx);

double mysecond();

static void
   free_and_null (char **ptr),
   usage         (char *progname);

int global_cnt=0;

pthread_mutex_t global_mtx1;
pthread_mutex_t global_mtx2;

FILE *logfile, *outputfile;

int
main (int argc, char *argv[])
{
   /* Declare and allocate space for the variables and arrays where we will
      store the optimization results including the status, objective value,
      and variable values. */

   int i,j,idx = 0, idxR = 0;                // Counters
   int ni, nj;  
   // 10x100        20x200      30x300  
   //  14.46952  19.329298   15.827

   double stepsize;
   double startM;
   //   double startM = -0.05;  //Try different value here

   int flag = 0;
   int status = 0;
   int nRows, nCols, nBins;
   double gap;
   char logbuffer[50];
   char data[100];
   double bestobjval, tempobjval;
   int nodecnt;
   int numfstbest;
   int mipitercnt = 0;
   double t0, t1;

   FILE *inputfile;
   inputfile = fopen(argv[1], "r");
   outputfile = fopen("../output/result.txt", "a");
   status = fscanf (inputfile," %d %d %lf %lf",  &ni, &nj, &startM, &stepsize);
   status = fscanf (inputfile, "%s", data);

   printf("%d %d %f %f\n", ni, nj, startM, stepsize);
   printf("%s\n", data);
   //   printf("Initial value of M is set to be %f\n", startM);

   sprintf(logbuffer, "../output/LogFile%dx%d_%3.2f_%3.2f.txt", ni, nj, startM, stepsize);
   // Parameters passed to the callback functions
   MINFO myMinfo;
   int **bigmRowIdx = NULL;      // row indices of the big Ms
   int **bigmColIdx = NULL;      // colume indices           
   double **bigmVals = NULL;

   double myBestObjVal = 1000000;

   // Solution parameters
   double objval;
   double   *x     = NULL;
   CPXENVptr     env = NULL;
   CPXLPptr      lp = NULL;
   double* lb = NULL;
   double* ub = NULL;

   // CPXFILEptr logfile = NULL; // create pointer for logfile

   /* Initialize the CPLEX environment */
   env = CPXopenCPLEX (&status);
   if ( env == NULL ) {
      char  errmsg[1024];
      fprintf (stderr, "Could not open CPLEX environment.\n");
      CPXgeterrorstring (env, status, errmsg);
      fprintf (stderr, "%s", errmsg);
      goto TERMINATE;
   }

   /* Setting mip options */
   status = DZturnoffpreproc(env);  if ( status )  goto TERMINATE;
   status = DZturnoffcuts(env);     if ( status )  goto TERMINATE;
   status = DZsetmipstrategy(env);  if ( status )  goto TERMINATE;
   status = DZsetmiplimits(env);    if ( status )  goto TERMINATE;
   status = DZsetmipoutput(env, logbuffer);    if ( status )  goto TERMINATE;
   //   fp = CPXfopen ("mylog.log", "w");

   CPXsetintparam(env,CPX_PARAM_PARALLELMODE, CPX_PARALLEL_DETERMINISTIC);
   if ( status )  goto TERMINATE;
   CPXsetintparam(env,CPX_PARAM_THREADS, 24);
   if ( status )  goto TERMINATE;

   /* Create the problem, using the filename as the problem name */
   lp = CPXcreateprob (env, &status, data);
   if ( lp == NULL ) {
      fprintf (stderr, "Failed to create LP.\n");
      goto TERMINATE;
   }

   status = CPXreadcopyprob (env, lp, data, NULL);
   if ( status ) {
      fprintf (stderr, "Failed to read and copy the problem data.\n");
      goto TERMINATE;
   }

   status = DZgetmipprobsize(env, lp, &nRows, &nCols, &nBins);
   if ( status )  goto TERMINATE;

   printf("nRows: %d  nCols: %d   nBins: %d \n", nRows, nCols, nBins);

   /* Allocate space for big M indices */
   bigmRowIdx = (int **) malloc(ni*sizeof(int *));
   bigmColIdx = (int **) malloc(ni*sizeof(int *));
   bigmVals =   (double **) malloc(ni*sizeof(double *));
   for(i = 0; i< ni; i++){
     bigmRowIdx[i] = (int *)malloc(nj*sizeof(int));
     bigmColIdx[i] = (int *)malloc(nj*sizeof(int));
     bigmVals[i]   = (double *)malloc(nj*sizeof(double));
   }

   /* Set the initial values of big M */
   for(i = 0; i< ni; i++){
     for(j=0; j<nj; j++)
       bigmVals[i][j] = startM;
   }

   /* Get the indices of big Ms in the A matrix, store 
      and check them */
   status = DZgetbigmindex(env, lp, ni, nj, bigmRowIdx, bigmColIdx);
   if ( status )  goto TERMINATE;

   /* Initialize cbhandle structure */
   myMinfo.flagPtr = &flag;
   myMinfo.ni = ni;
   myMinfo.nj = nj;
   myMinfo.bigmRowIdx = bigmRowIdx;
   myMinfo.bigmColIdx = bigmColIdx;
   myMinfo.myBestObjValPtr = &myBestObjVal;
   myMinfo.bigmValsPtr     = &bigmVals;
   myMinfo.stepsize = stepsize;


   /* Initialize the mutex */
   pthread_mutex_init(&myMinfo.mtx,NULL);
   pthread_mutex_init(&myMinfo.mtx1,NULL);
   pthread_mutex_init(&global_mtx1,NULL);
   pthread_mutex_init(&global_mtx2,NULL);

   /* Set up to use MIP solve callbacks */
   status = CPXsetsolvecallbackfunc (env, usersetM, &myMinfo);
   if ( status ) {
     fprintf (stderr,"Failed to install sovle callback functions.\n");
     goto TERMINATE;
   }

   /* Set up to use MIP info callbacks */
   status = CPXsetmipcallbackfunc (env, flagcallback, &myMinfo);
   if ( status ) {
     fprintf (stderr,"Failed to install info callback functions.\n");
     goto TERMINATE;
   }


   t0 = mysecond();

   /* Checking , changing coefficients and resolve */
   do {
     flag = 0; //reset the flag
     mipitercnt++;
     //printf("This is another iteration\n");

     for(i=0; i<ni; i++){
       for(j=0; j<nj; j++)
	 status = CPXchgcoef(env, lp, bigmRowIdx[i][j], bigmColIdx[i][j], bigmVals[i][j]);
     }

     /* Optimize the problem and obtain the solution */
     status = CPXmipopt (env, lp);
     if ( status ) {
       fprintf (stderr,"Failed to optimize MIP.\n");
       goto TERMINATE;
     }

     status  = CPXgetmipobjval (env, lp, &objval);

     //Store the numfstbest
     if(objval != tempobjval){
       tempobjval = objval;
       numfstbest = mipitercnt;
     }

     status  = CPXgetbestobjval (env, lp, &bestobjval);
     //     printf("Best upper bound: %12.6g   Best lower bound: %12.6g\n", objval, bestobjval);
   } while(flag == 1);

   t1 = mysecond();

   /* Get the optimal cost */
   status  = CPXgetmipobjval (env, lp, &objval);
   if ( status ) {
     fprintf (stderr,"Failed to obtain objective value.\n");
     goto TERMINATE;
   }
   printf("The objective value is: %f\n", objval);

   nodecnt = CPXgetnodecnt (env, lp);
   printf("The number of nodes is %d\n", nodecnt);

   status  = CPXgetbestobjval (env, lp, &bestobjval);
   if ( status ) {
     fprintf (stderr,"Failed to obtain objective value.\n");
     goto TERMINATE;
   }
   printf("The best objective value is: %f\n", bestobjval);
   printf("Number of iteration is: %d\n", mipitercnt);
   printf("Elaspsed time is: %f\n", t1-t0);

   /* Allocate space for solution */
   x = (double *) malloc (nCols*sizeof(double));
   if ( x == NULL ) {
     fprintf (stderr, "No memory for solution values.\n");
     goto TERMINATE;
   }

   /* Get the optimal solution */
   status = CPXgetmipx (env, lp, x, 0, nCols-1);
   if ( status ) {
     fprintf (stderr, "Failed to obtain solution.\n");
     goto TERMINATE;
   }

   //Here is the relationship between nCols and ni nj
   // nCols = (ni+1)*nj + ni + 1
   // The last ni is values of Y
   // There are (ni+1)*nj Xij and first variable is Z
   for (i=(ni+1)*nj+1; i<(ni+1)*nj + ni+1; i++){
     printf("%12.6f \n",x[i]) ;
   }

   fprintf(logfile,"X values\n");
   for(i= 0; i<ni; i++){
     for(j=0; j<nj; j++){
       idx = (ni + nj*i + (j+1)) - ni; // index of x at the solution vector
       idxR = (ni+1)*nj + i+1;         // index of y(b) at the solution vector
       gap = x[idx] + bigmVals[i][j]*x[idxR];

       if(x[idxR] <= 1e-8)
	 fprintf(logfile,"%12.6f %12.6f %12.6f      ------\n", x[idx], bigmVals[i][j], x[idxR]);
       else
	 fprintf(logfile,"%12.6f %12.6f %12.6f %12.6f\n", x[idx], bigmVals[i][j], x[idxR], gap);
     }
   }

   fprintf(logfile,"Y values\n");
   for(i= nCols-ni; i<nCols; i++){
     fprintf(logfile,"%12.6f \n", x[i]);
   }

   fprintf(logfile,"MIP iterations counts: %d\n", mipitercnt);

   /* Set near-zero values to 0 */
   for (int j = 0; j < nCols; j++) {
   	  if( x[j] > 0 && x[j] < 0.000000001 ) x[j] = 0.0;
   }

   /* Print out formatted result to the output file */
   fprintf(outputfile, "%6d  %6d  %6.4f %6.4f  ",ni, nj, startM, stepsize);   //echo ni nj startM delta
   fprintf(outputfile, "%6d  %6d  %6d  %12.6f %12.6f\n", numfstbest, mipitercnt, nodecnt, t1-t0, objval);
   
TERMINATE:

   /* Free up the solution */
   free_and_null ((char **) &x);
   free_and_null ((char **) &bigmRowIdx);
   free_and_null ((char **) &bigmColIdx);
   free_and_null ((char **) &bigmVals);
   pthread_mutex_destroy (&myMinfo.mtx);
   pthread_mutex_destroy (&myMinfo.mtx1);
   pthread_mutex_destroy (&global_mtx1);
   pthread_mutex_destroy (&global_mtx2);

   /* Free up the problem as allocated by CPXcreateprob, if necessary */
   if ( lp != NULL ) {
      status = CPXfreeprob (env, &lp);
      if ( status ) {
         fprintf (stderr, "CPXfreeprob failed, error code %d.\n", status);
      }
   }

   /* Free up the CPLEX environment, if necessary */
   if ( env != NULL ) {
      status = CPXcloseCPLEX (&env);
      if ( status ) {
      char  errmsg[1024];
         fprintf (stderr, "Could not close CPLEX environment.\n");
         CPXgeterrorstring (env, status, errmsg);
         fprintf (stderr, "%s", errmsg);
      }
   }
   
   if ( logfile != NULL )  CPXfclose (logfile); // close logfile
   fclose(inputfile);
   fclose(outputfile);
   return (status);

}  /* END main */


static int CPXPUBLIC
usersetM      (CPXCENVptr   env,
               void         *cbdata,
               int          wherefrom,
               void         *cbhandle,
               int          *useraction_p
	       )   
{
   int status = 0;
   int nodecount,seqnum, itercnt=0;
   int rows, cols;
   int i,j,idx,idxR, idxM;
   int mvalidind = 1;
   double coef, objval, bestobjval, bestintval,cutoffval;
   double tol = 0.0000000001;
   char lu[1];
   double bd[1];
   int id[1];
   int depth;
   double   *x   = NULL;
   char nodeLPName[50];
   //   int i_min = 0;
   //   int i_max = 30;
   int threadnum, numthreads;
   int cbindIdct = 0;
   int rbindIdct = 0;
   int ybindIdct = 0;
   int      *cstat = NULL;
   int      *rstat = NULL;
   double mingap = -5;
   global_cnt += 1; 
   //printf("%d %d\n", pthread_self(), global_cnt);
 
 
   void *status_p;

   MINFOptr info = (MINFOptr) cbhandle;   //cast the cbhandle

   /* Read in info of coordinates of big M */
   int ni = info -> ni;                    // Read 
   int nj = info -> nj;                    // Read
   double stepsize = info->stepsize;       // Read
   int **bigmRowIdx = info -> bigmRowIdx;  // Read
   int **bigmColIdx = info -> bigmColIdx;  // Read

   int *flagPtr = info->flagPtr;                       // Global Write
   double *myBestObjValPtr = info -> myBestObjValPtr;  // Global Write
   double ***bigmValsPtr = info -> bigmValsPtr;        // Global Write

   CPXLPptr nodelp;      //Can be modified
   CPXCLPptr lp;         //Constant pointer. Can not be modified
   
   //   *useraction_p = CPX_CALLBACK_DEFAULT;

   /* Get pointer to root LP problems */
   status = CPXgetcallbacklp(env, cbdata, wherefrom, &lp);
   if ( status )  goto TERMINATE;

   /* Get pointer to node LP problems */
   status = CPXgetcallbacknodelp (env, cbdata, wherefrom, &nodelp);
   if ( status )  goto TERMINATE;
   //   printf("%d\n", lp);       
   //   printf("%d\n", nodelp);

   /* Find out what nodes is being processed */
   status = CPXgetcallbackinfo (env, cbdata, wherefrom, CPX_CALLBACK_INFO_NODE_COUNT, &nodecount);
   if ( status )  goto TERMINATE;

   status = CPXgetcallbackinfo (env, cbdata, wherefrom, CPX_CALLBACK_INFO_MY_THREAD_NUM, &threadnum);
   if ( status )  goto TERMINATE;

   status = CPXgetcallbackinfo (env, cbdata, wherefrom, CPX_CALLBACK_INFO_USER_THREADS, &numthreads);
   if ( status )  goto TERMINATE;
   


   /* Find out the node sequential number */
   //   status = CPXgetcallbacknodeinfo( env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_SEQNUM, &seqnum);
   if ( status )  goto TERMINATE;
   //   status = CPXgetcallbacknodeinfo( env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_DEPTH, &depth);
   if ( status )  goto TERMINATE;
   //   printf("I am at node %d of depth %d \n", seqnum, depth);

   /*
   if(seqnum <= i_max && seqnum > i_min){
     sprintf(nodeLPName, "tempnodelp%d.lp",seqnum);
     printf("%s\n", nodeLPName);
     CPXwriteprob(env, nodelp, nodeLPName,NULL);
   }
   */

   //   if (nodecount < 1) goto TERMINATE; //skip the solving process
   //   if(nodecount == 1)
   //     printf("I am going off the root node!\n");

   //   cols = CPXgetnumcols (env, nodelp);

   cols = CPXgetnumcols (env, lp);
   if ( cols <= 0 ) {
     fprintf (stdout, "Can't get number of columns.\n");
     status = CPXERR_CALLBACK;
     goto TERMINATE;
   }

   rows = CPXgetnumrows (env, lp);
   if ( rows <= 0 ) {
     fprintf (stdout, "Can't get number of rows.\n");
     status = CPXERR_CALLBACK;
     goto TERMINATE;
   }

   x = (double *) malloc (cols * sizeof (double));
   cstat = (int *) malloc (cols*sizeof(int));
   rstat = (int *) malloc (rows*sizeof(int));

   /* Solve initial node with barrier. others with dual */
   //   status = CPXprimopt(env,nodelp);
   //status = CPXdualopt(env,nodelp);
   if (nodecount < 1)   status = CPXprimopt(env, nodelp);
   else                 status = CPXdualopt(env, nodelp);
   //     status = CPXprimopt(env, nodelp);

   if ( status ) {
     fprintf (stdout, "Can't solve the node LP.");
     goto TERMINATE;
   }

   /* Retrieve the solution at the node LP */
   status = CPXgetx(env, nodelp, x, 0, cols-1);
   if ( status ) {
     fprintf (stdout, "Can't get node solution.");
     goto TERMINATE;
   }

   //      rstat = (int *) malloc (*sizeof(int));
   if ( cstat == NULL ) {
     fprintf (stderr, "No memory for basis statuses.\n");
     goto TERMINATE;
   }
   //      printf("%d\n", cstat);

   status = CPXgetbase (env,  nodelp, cstat, rstat);
   //   printf("%d\n", cstat[5]);
   if ( status ) {
     fprintf (stderr, "Failed to get basis; error %d.\n", status);
     goto TERMINATE;
   }

   status = CPXgetobjval(env, nodelp, &objval);
   if ( status ) {
     fprintf (stdout, "Can't get objective value of node LP.");
     goto TERMINATE;
   }

   // Callback best objective value among the remaing nodes
   //   printf("I am good here!\n");
   status = CPXgetcallbackinfo (env, cbdata, wherefrom,
                                   CPX_CALLBACK_INFO_BEST_REMAINING, &bestobjval);
   if ( status )  goto TERMINATE;
   status = CPXgetcallbackinfo (env, cbdata, wherefrom,
                                   CPX_CALLBACK_INFO_BEST_INTEGER, &bestintval);
   if ( status )  goto TERMINATE;
   status = CPXgetcallbackinfo (env, cbdata, wherefrom,
                                   CPX_CALLBACK_INFO_CUTOFF, &cutoffval);
   if ( status )  goto TERMINATE;


   if(objval >0 && objval < *myBestObjValPtr)
     *myBestObjValPtr = objval;
   //   printf("I am at thread %d node %d.\n", threadnum, seqnum);

   // Apply mutex to change the coefficient matrix
   pthread_mutex_lock (&info->mtx);
   for(i = 0; i< ni; i++){
     for(j = 0; j< nj; j++){

       idx = (ni + nj*i + (j+1)) - ni; // index of x at the solution vector    /why this is correct?
       idxR = (ni+1)*nj + i+1;         // index of y(b) at the solution vector /why this is correct?

       status = CPXgetcoef(env, lp, bigmRowIdx[i][j], bigmColIdx[i][j], &coef);

       if(status == 1){
	 printf("I have something wrong here!\n");
	 goto TERMINATE;
       }
       
       rbindIdct = rstat[bigmRowIdx[i][j]];
       cbindIdct = cstat[idx];
       ybindIdct = cstat[idxR];

       // Big M Checking routine
       // Binding
      if( coef + 1.0 <= tol || cstat[idx] == CPX_FREE_SUPER ||rbindIdct == CPX_BASIC) {
	 itercnt ++; 
	 continue;
       }

      // binding then increase the big-M
       else{
	 if(coef - stepsize >= -1.0) coef -= stepsize;
	 else coef = -1.0;

	 (*bigmValsPtr)[i][j] = coef;
       }
     }
   }

  // End of checking loop
   pthread_mutex_unlock (&info->mtx);

   //   printf("itercnt is %d status is %d\n", itercnt, status);
  if (itercnt == ni*nj){
     //    printf("mip passed at thread %d.\n", threadnum);
     goto TERMINATE;
   }
  else{//Go to the next iteration
    //    printf(" Go to the next iteration! \n");
    //fprintf(logfile, "Best lower bound is %12.6f \n", cutoffval);
    //    printf("mip failed at thread %d of %d.\n", threadnum, numthreads);
    //    free_and_null ((char **) &x);

    pthread_mutex_lock(&info->mtx1);
    *flagPtr = 1;
    pthread_mutex_unlock (&info->mtx1);
    printf("mip failed at id %d \n", pthread_self());
    //    printf("flag is %d \n", *flagPtr);
    goto TERMINATE;
  }

 TERMINATE:
   *useraction_p = CPX_CALLBACK_SET;
   free_and_null ((char **) &x);
   free_and_null ((char **) &rstat);
   free_and_null ((char **) &cstat);
   return (status);

}   //End of callback function usersetM



static int flagcallback (CPXCENVptr env, void *cbdata, int wherefrom,
                   void *cbhandle)
{
  int status = 0;
  int threadnum;

  //  printf(" I am here !\n");
  MINFOptr info = (MINFOptr) cbhandle;
  int *flagPtr = info->flagPtr;                       // Global Write

  CPXgetcallbackinfo (env, cbdata, wherefrom, CPX_CALLBACK_INFO_MY_THREAD_NUM, &threadnum);
  //  printf(" I am at thread %d\n", threadnum);
  //  printf("flag here is %d.\n", *flagPtr);
 
  if(*flagPtr == 1){
    printf("%d\n", pthread_self());
    status = 1;
  }
  return status;
} 

/* This simple routine frees up the pointer *ptr, and sets *ptr to NULL */

static void
free_and_null (char **ptr)
{
   if ( *ptr != NULL ) {
      free (*ptr);
      *ptr = NULL;
   }
} /* END free_and_null */ 


static void
usage (char *progname)
{
   fprintf (stderr,"Usage: %s filename\n", progname);
   fprintf (stderr,"   where filename is a file with extension \n");
   fprintf (stderr,"      MPS, SAV, or LP (lower case is allowed)\n");
   fprintf (stderr,"  This program uses the CPLEX MIP optimizer.\n");
   fprintf (stderr," Exiting...\n");
} /* END usage */


int DZturnoffpreproc(CPXENVptr env){
   int status = 0;

  /* Turn on or off preprocessor and presolve for the nodes */
   status = CPXsetintparam(env, CPX_PARAM_PREIND, CPX_OFF);       // presolve 
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_PRESLVND, CPX_OFF);     // presolve at the node
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_RELAXPREIND, CPX_OFF);  // presolve of relaxed lp at the root node  
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_COEREDIND, CPX_OFF);   // coefficient reduction
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_AGGIND, CPX_OFF);      // aggregator
   if ( status )  return status;

   /* Assure linear mappings between the presolved and original
      models */
   status = CPXsetintparam (env, CPX_PARAM_PRELINEAR, 0);
   if ( status )  return status;

   /* Let MIP callbacks work on the original model */
   CPXsetintparam(env, CPX_PARAM_MIPCBREDLP, CPX_OFF);
   if ( status )  return status;

   return 0;
}

int DZturnoffcuts(CPXENVptr env){
  int status = 0;
   /* Set cut options */
   status = CPXsetintparam(env,CPX_PARAM_CLIQUES , -1);	    if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_COVERS , -1);	    if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_DISJCUTS , -1);	    if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_FLOWPATHS , -1);    if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_GUBCOVERS , -1);    if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_FRACCUTS , -1);	    if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_IMPLBD , -1);       if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_FLOWCOVERS , -1);   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_MIRCUTS , -1 );    if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_ZEROHALFCUTS ,-1 );   if ( status ) return status;
   return 0;
}

int DZsetmipstrategy(CPXENVptr env){

  int status;
   /* Select the root algorithm */
   status = CPXsetintparam (env, CPX_PARAM_STARTALG,CPX_ALG_AUTOMATIC);
   if ( status )  return status;

   /* Set controlling strategies */
   status = CPXsetintparam(env, CPX_PARAM_MIPEMPHASIS, 0);	 //3 to emphasize the best bound
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_NODESEL, 1); //or 2
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_VARSEL, 0); //0[Auto]
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_BRDIR, 0); //-1[Downward] 0[Auto] 1[Upward]
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_PROBE, -1);
   if ( status )  return status;
   status = CPXsetintparam(env, CPX_PARAM_ADVIND, 1);  //Need to be turned off otherwise it will cause trouble for the next iteration
   if ( status )  return status;
   /* Set heuristic strategy */
   status = CPXsetintparam(env, CPX_PARAM_HEURFREQ , -1 );    //-1 to turn it off
   if ( status )  return status;
   //   status = CPXsetintparam(env, CPX_PARAM_BBINTERVAL , 0);   
   //   if ( status )  return status;
   //   CPXsetintparam(env, CPX_PARAM_FPHEUR , -1 );    
   //   CPXsetintparam(env, CPX_PARAM_LBHEUR , 1 );    

   /* Turn on traditional search for user with control callbacks */
   //   status = CPXsetintparam (env, CPX_PARAM_MIPSEARCH, CPX_MIPSEARCH_TRADITIONAL);
   //   if ( status )  return status;

   /* Set parallel mode and threads */
   /*
   CPXsetintparam(env,CPX_PARAM_PARALLELMODE, CPX_PARALLEL_DETERMINISTIC);
   if ( status )  return status;
   CPXsetintparam(env,CPX_PARAM_THREADS, 8);
   if ( status )  return status;
   */
   return 0;

}

int  DZsetmiplimits  (CPXENVptr env){
  int status = 0;
/* Set terminating criterion*/
     CPXsetdblparam(env, CPX_PARAM_EPGAP, 0.00000 );   
     CPXsetdblparam(env, CPX_PARAM_EPAGAP, 0.00000 );   
     CPXsetintparam(env, CPX_PARAM_INTSOLLIM, 1000000000); 
  //   if ( status )  return status;
  //   CPXsetintparam(env, CPX_PARAM_NODELIM, 10000000); 
  //   if ( status )  return status;
   //   CPXsetdblparam(env,  CPX_PARAM_TILIM, 3600);    
   CPXsetintparam(env,  CPX_PARAM_CLOCKTYPE, 1); 
   // CPXsetintparam(env,CPX_PARAM_MIPSTART, 0); 
   return status;
}


int DZsetmipoutput (CPXENVptr env, char *logbuffer){
  int status = 0;

  /* Turn on output to the screen */
   status = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON);
   if ( status ) {
      fprintf (stderr, 
               "Failure to turn on screen indicator, error %d.\n", status);
      return status;
   }

   /* Set output options */
   status = CPXsetintparam(env,CPX_PARAM_MIPINTERVAL , 50);
   if ( status )  return status;
   status = CPXsetintparam(env,CPX_PARAM_MIPDISPLAY , 2);
   if ( status )  return status;


   // ======= Direct output to a logfile
   logfile = CPXfopen (logbuffer, "w");
   if ( logfile == NULL )  return status;
   status = CPXsetlogfile (env, logfile);
   if ( status )  return status;

   fprintf(logfile,"\n-------------\n");

   return status;
}



int DZgetmipprobsize(CPXENVptr env,CPXLPptr lp, int *nRowsPtr, int *nColsPtr, int *nBinsPtr){

  int status = 0;
  *nRowsPtr = CPXgetnumrows (env, lp);
  *nColsPtr = CPXgetnumcols (env, lp);
  *nBinsPtr = CPXgetnumbin(env, lp);
  //   printf("Num rows and colums: %d %d \n",nRows, nCols);

  //   printf("Number of binary variables are: %d \n", nBins);

  return status;
}

int DZgetbigmindex(CPXENVptr env, CPXLPptr lp, int ni, int nj, 
		  int **bigmRowIdx, int **bigmColIdx) {
  int status = 0;
  int i,j,k;
  char tempColName[10];
  char tempRowName[10];
  double coef;

  for(i = 0; i<ni; i++){
     for(j = 0; j<nj; j++){
       sprintf(tempRowName, "c%d", i*nj+j+2);
       sprintf(tempColName, "b%d", i+1);

       CPXgetrowindex(env,lp,tempRowName, &bigmRowIdx[i][j]);
       CPXgetcolindex(env,lp,tempColName, &bigmColIdx[i][j]);
   
       status = CPXgetcoef(env, lp, bigmRowIdx[i][j], bigmColIdx[i][j], &coef);
       if(coef != -1){
	 fprintf (stderr, "Big M checking faied.\n");
	 return 1;
       }
     }
   }
   printf("All the coordinates of big M are checked!\n");
   return status;
}

double mysecond()
{
   struct timeval tp;
   struct timezone tzp;
   int i;
   i = gettimeofday(&tp,&tzp);
   return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
}

